Graphical Analysis

1 📖 Main Menu

using StatsBase, DataFrames, DataFrameMacros
using RCall, CategoricalArrays, TexTables
using Distributions, Gadfly, Compose
using CSV, RDatasets, MLJ
using MLJ: schema

include("pubh.jl")
@rlibrary readr
@rimport pubh

2 Data

The Western Collaborative Group Study (WCGS) is a well known prospective cohort study. Male participants aged 39 to 59 from 10 California companies were originally selected to study the relationship between behaviour pattern and the risk of coronary heart disease (CHD).

wcgs = read_rds("data/wcgs.rds") |> rcopy
wcgs |> schema

┌──────────┬───────────────────────────────┬────────────────────────────────────
│ names    │ scitypes                      │ types                             ⋯
├──────────┼───────────────────────────────┼────────────────────────────────────
│ id       │ Count                         │ Int64                             ⋯
│ age      │ Count                         │ Int64                             ⋯
│ height   │ Continuous                    │ Float64                           ⋯
│ weight   │ Continuous                    │ Float64                           ⋯
│ sbp      │ Count                         │ Int64                             ⋯
│ dbp      │ Count                         │ Int64                             ⋯
│ chol     │ Union{Missing, Count}         │ Union{Missing, Int64}             ⋯
│ beh_pat  │ Multiclass{4}                 │ CategoricalValue{String, UInt32}  ⋯
│ ncigs    │ Count                         │ Int64                             ⋯
│ dib_pat  │ Multiclass{2}                 │ CategoricalValue{String, UInt32}  ⋯
│ chd      │ Multiclass{2}                 │ CategoricalValue{String, UInt32}  ⋯
│ type_chd │ Multiclass{4}                 │ CategoricalValue{String, UInt32}  ⋯
│ time     │ Count                         │ Int64                             ⋯
│ arcus    │ Union{Missing, Multiclass{2}} │ Union{Missing, CategoricalValue{S ⋯
│ smoker   │ Multiclass{2}                 │ CategoricalValue{String, UInt32}  ⋯
└──────────┴───────────────────────────────┴────────────────────────────────────
                                                                1 column omitted

The kfm data frame was collected by Kim Fleischer Michaelsen and contains data for 50 infants of age approximately 2 months. They were weighed immediately before and after each breast feeding. and the measured intake of breast milk was registered along with various other data.

kfm = read_rds("data/kfm.rds") |> rcopy
kfm |> schema

┌────────────┬───────────────┬──────────────────────────────────┐
│ names      │ scitypes      │ types                            │
├────────────┼───────────────┼──────────────────────────────────┤
│ no         │ Count         │ Int64                            │
│ dl_milk    │ Continuous    │ Float64                          │
│ sex        │ Multiclass{2} │ CategoricalValue{String, UInt32} │
│ weight     │ Continuous    │ Float64                          │
│ ml_suppl   │ Count         │ Int64                            │
│ mat_weight │ Count         │ Int64                            │
│ mat_height │ Count         │ Int64                            │
└────────────┴───────────────┴──────────────────────────────────┘

Data on reported cases of influenza by age group in the 1957 pandemic in England and Wales.

flu = DataFrame(CSV.File("data/fluraw.csv"))
flu |> schema

┌───────┬────────────────┬───────┐
│ names │ scitypes       │ types │
├───────┼────────────────┼───────┤
│ week  │ ScientificDate │ Date  │
│ child │ Count          │ Int64 │
│ young │ Count          │ Int64 │
│ mid   │ Count          │ Int64 │
│ old   │ Count          │ Int64 │
└───────┴────────────────┴───────┘

Data on birth weights from newborns in a London Hospital.

birth = read_rds("data/birthwt.rds") |> rcopy
birth |> schema

┌───────┬───────────────┬──────────────────────────────────┐
│ names │ scitypes      │ types                            │
├───────┼───────────────┼──────────────────────────────────┤
│ low   │ Multiclass{2} │ CategoricalValue{String, UInt32} │
│ age   │ Count         │ Int64                            │
│ lwt   │ Continuous    │ Float64                          │
│ race  │ Multiclass{3} │ CategoricalValue{String, UInt32} │
│ smoke │ Multiclass{2} │ CategoricalValue{String, UInt32} │
│ ptl   │ Count         │ Int64                            │
│ ht    │ Multiclass{2} │ CategoricalValue{String, UInt32} │
│ ui    │ Multiclass{2} │ CategoricalValue{String, UInt32} │
│ ftv   │ Count         │ Int64                            │
│ bwt   │ Count         │ Int64                            │
│ Race  │ Multiclass{2} │ CategoricalValue{String, UInt32} │
└───────┴───────────────┴──────────────────────────────────┘

Data on the energy expenditure in groups of lean and obese women.

energy = rcopy(R"ISwR::energy")
energy |> schema

┌─────────┬───────────────┬──────────────────────────────────┐
│ names   │ scitypes      │ types                            │
├─────────┼───────────────┼──────────────────────────────────┤
│ expend  │ Continuous    │ Float64                          │
│ stature │ Multiclass{2} │ CategoricalValue{String, UInt32} │
└─────────┴───────────────┴──────────────────────────────────┘

Fisher’s iris dataset on measurements for three species of iris.

iris = rcopy(R"datasets::iris")
iris |> schema

┌──────────────┬───────────────┬──────────────────────────────────┐
│ names        │ scitypes      │ types                            │
├──────────────┼───────────────┼──────────────────────────────────┤
│ Sepal_Length │ Continuous    │ Float64                          │
│ Sepal_Width  │ Continuous    │ Float64                          │
│ Petal_Length │ Continuous    │ Float64                          │
│ Petal_Width  │ Continuous    │ Float64                          │
│ Species      │ Multiclass{3} │ CategoricalValue{String, UInt32} │
└──────────────┴───────────────┴──────────────────────────────────┘

3 Distributions

3.1 Histograms

We use histograms to look at the distribution of continuous variables. Histograms make most sense for relatively large data (\(\geq 100\) observations). In this fist histogram, we show the density.

plot(
  birth,
  x = :bwt, Guide.xlabel("Birth weight (g)"),
  Guide.ylabel("Density"),
  Geom.histogram(bincount=30, density=true)
)

Figure 1: Distribution of birth weight from newborns in a London Hospital.

The default for histograms is to show the absolute frequency (counts):

plot(
  wcgs,
  x = :sbp, 
  Geom.histogram(bincount=30),
  Guide.xlabel("SBP (mm Hg)"),
  Guide.ylabel("Counts"),
  Coord.cartesian(xmin=80, ymax=700)
)

Figure 2: Distribution of systolic blood pressure (SBP) from subjects in the WCGS dataset.

We can compare distributions using the density or the probability density function (pdf), as shown here:

plot(
  wcgs,
  x = :sbp, xgroup = :chd,
  Geom.subplot_grid(
    Geom.histogram(bincount=30, density=true)
    ),
  Guide.xlabel("SBP (mm Hg)"),
  Guide.ylabel("Density")
)

Figure 3: Comparison of systolic blood pressure (SBP) distributions by coronary heart disease (CHD) events.

3.2 Density Plots

Another way to look at the distribution of a continuous variable is with density plots.

plot(
  wcgs,
  x = :sbp, color = :chd,
  Geom.density,
  Guide.xlabel("SBP (mm Hg)"),
  Guide.ylabel("Density"),
  Guide.colorkey(title="CHD"),
  style(line_width=1mm),
  Coord.cartesian(xmin=80)
)

Figure 4: Distribution of systolic blood pressure (SBP) by coronary heart disease (CHD) event.

Diastolic blood pressure (DBP) is skewed to the right (see the QQ-Plot). In the following example, we log-transform DBP.

plot(
  wcgs,
  x = log.(wcgs.dbp), color = :chd,
  Geom.density,
  Guide.xlabel("log (DBP)"),
  Guide.ylabel("Density"),
  Guide.colorkey(title="CHD"),
  style(line_width=1mm),
)

Figure 5: Comparison of diastolic blood pressure (DBP) distributions between coronary heart disease (CHD) status of the subjects. DBP was log-transformed to bring it closer to a normal distribution.

Tip

We have missing values in chol, thus we need to drop them before constructing a density plot or else, we would have an error.

plot(
  @subset(wcgs, !ismissing(:chol)),
  x = :chol, color = :chd,
  Geom.density,
  Guide.xlabel("Cholesterol (mg/dl)"),
  Guide.ylabel("Density"),
  Guide.colorkey(title="CHD"),
  style(line_width=1mm),
  Coord.cartesian(ymax=0.011)
)

Figure 6: Distribution of blood cholesterol by coronary heart disease (CHD) event from the WCGS dataset.

3.3 QQ-Plots

The best way to determine if a continuous variable is normally distributed or not is with quantile-quantile plots (QQ-plots). We plot the quantiles of our variable of interest against quantiles from the standard normal distribution (which has a mean \(\mu=0\) and a standard deviation \(\sigma=1\)). This type of QQ-plots against the normal distribution are known as QQ-normal plots. If the variable is normally distributed, then a linear relationship will be observed.

qq_plot(birth.bwt; ylab = "Birth weight (g)")

Figure 7: QQ-plot of birth weights from newborns in a London hospital.

Note

The function qq_plot is part of pubh.jl, a wrapper to functions from Gadfly constructs QQ-plots.

Right-skewed variables, show an upper right curve in the QQ-plot:

qq_plot(wcgs.dbp; ylab="DBP (mm Hg)")

Figure 8: QQ-plot of diastolic blood pressure (DBP) in subjects from the WCGS dataset.

The next plot, shows the use of log scales.

plot(
  wcgs,
  x=fit.([Normal], [wcgs.dbp]), y=:dbp,
  Stat.qq,
  Guide.xlabel("Theoretical quantiles"),
  Guide.ylabel("DBP (mm Hg)"),
  Scale.y_log10,
  Theme(default_color="IndianRed")
)

Figure 9: QQ-plot of diastolic blood pressure (DBP) in subjects from the WCGS dataset.

We can plot the log of the variable, instead of using log-scales.

qq_plot(log.(wcgs.dbp), ylab="log (DBP)")

4 Associations Between Continuous Variables

4.1 Scatter Plots

We use scatter plots to look at the relationship between two continuous variables. By default, the dependent variable (response) is plotted on the \(y\)-axis while the independent variable (explanatory) is plotted on the \(x\)-axis.

From the kfm dataset, let’s see if there is a relationship between the weight of the mother and the breast-milk intake of the child.

plot(
  kfm,
  x = :mat_weight, y = :dl_milk, color = :sex,
  Guide.xlabel("Maternal weight (kg)"),
  Guide.ylabel("Breast milk intake (dl/day) ")
)

Figure 10: Relationship between maternal weight and breast-milk intake in infants.

Note

We can extend plots by adding new layers. In the following example, we use two layers: one for the scatter plot and a second to add a line representing a linear fit.

plot(
  kfm,
  x = :mat_weight, y = :dl_milk, color = :sex,
  Geom.point, 
  layer(
    Stat.smooth(method=:lm), 
    Geom.line, 
    Geom.ribbon
    ),
    Guide.xlabel("Maternal weight (kg)"),
    Guide.ylabel("Breast milk intake (dl/day) ")
)

Figure 11: Relationship between maternal weight and breast-milk intake in infants.

In the next plot, instead of showing the linear fit, we add a loess trend with a smoother:

plot(
  kfm,
  x = :mat_weight, y = :dl_milk, xgroup = :sex,
  Geom.subplot_grid(Geom.point, Geom.smooth(smoothing=0.9)),
  Guide.xlabel("Maternal weight (kg)"),
  Guide.ylabel("Breast milk intake (dl/day) ")
)

Figure 12: Relationship between maternal weight and breast-milk intake in infants.

In the next plot, we use different line colours and different symbols (by sex):

plot(
  kfm,
  x = :mat_weight, y = :dl_milk, 
  color = :sex, shape = :sex,
  Geom.point, Geom.smooth(method = :lm),
  Guide.xlabel("Maternal weight (kg)"),
  Guide.ylabel("Breast milk intake (dl/day) ")
)

Figure 13: Relationship between maternal weight and breast-milk intake in infants.

4.2 Line Charts

In this example, our dataset is in wide format.

Reshaping wide to long data:

flu_melt = stack(flu, Not(:week))
flu_melt |> head

5×3 DataFrame

Row	week	variable	value
	Date	String	Int64
1	1957-08-25	child	0
2	1957-09-01	child	0
3	1957-09-08	child	0
4	1957-09-15	child	23
5	1957-09-22	child	63

plot(
  flu_melt,
  x = :week, y = :value, color = :variable,
  Geom.point, Geom.path,
  Guide.xlabel("Date"),
  Guide.ylabel("Number of cases"),
  Guide.colorkey("Age group")
)

Figure 14: Dynamics of reported influenza cases by age group, in the 1957 pandemic in England and Wales.

5 Comparing Groups

5.1 Box-Plots

When we are comparing continuous variables, between two or more groups, box plots are the best option, particularly if the number of observations in the groups is relatively large \((\geq 30)\).

plot(
  @subset(wcgs, !ismissing(:chol)),
  x = :chd, y = :chol,
  Geom.boxplot,
  Guide.xlabel("Coronary Heart Disease"),
  Guide.ylabel("Cholesterol (mg/dl)")
)

Figure 15: Box plots of blood cholesterol by coronary heart disease (CHD) event from the WCGS dataset.

In the previous figure the presence of an outlier is clear. If we would like to remove that outlier, we would need to declare that in the report. For demonstration purposes, if we do not want to show the outlier in the plot, we have the option to filter the data using the @subset command.

plot(
  @subset(wcgs |> dropmissing, :chol .< 500),
  x = :chd, y = :chol,
  Geom.boxplot,
  Guide.xlabel("Coronary Heart Disease"),
  Guide.ylabel("Cholesterol (mg/dl)"),
  Theme(default_color="MidnightBlue")
)

plot(
  birth,
  x = :race, y = :bwt, color = :smoke,
  Geom.boxplot,
  Guide.xlabel(""),
  Guide.ylabel("Birth weight (g)"),
  Guide.colorkey("")
)

Figure 16: Comparison of birth weight distributions by ethnicity and smoking status of the mother.

Faceting:

plot(
  birth,
  x = :Race, y = :bwt, color=:Race, xgroup = :smoke,
  Geom.subplot_grid(Geom.boxplot),
  Guide.xlabel(""),
  Guide.ylabel("Birth weight (g)"),
  Guide.colorkey("")
)

Figure 17: Comparison of birth weight distributions by ethnicity and smoking status of the mother.

5.2 Strip Charts

When we want to compare groups and the number of observations is relatively small \((n<30)\), strip charts are superior than box-plots as we can show all observations.

strip_error(
  energy, "stature", "expend",
  xlab = "Stature",
  ylab = "Energy expenditure (MJ)"
)

Figure 18: Comparison of energy expenditure between obese and lean women.

5.3 Error Bars

birth_bst = combine(groupby(birth, [:smoke, :race]), :bwt => cis => AsTable)

6×5 DataFrame

Row	smoke	race	outcome	lower	upper
	Cat…	Cat…	Float64	Float64	Float64
1	Non-smoker	White	3428.75	3218.93	3638.57
2	Non-smoker	African-American	2854.5	2550.09	3158.91
3	Non-smoker	Other	2815.78	2628.31	3003.25
4	Smoker	White	2826.85	2656.57	2997.12
5	Smoker	African-American	2504.0	2109.15	2898.85
6	Smoker	Other	2757.17	2298.84	3215.49

plot(
  birth_bst,
  x=:race, y=:outcome, xgroup=:smoke,
  ymin=:lower, ymax=:upper,
  Geom.subplot_grid(
    Geom.point,
    Geom.errorbar
    ),
  Guide.ylabel("Birth weight (g)"),
  Guide.xlabel(""),
  Theme(default_color="Firebrick")
)

Figure 19: Comparison of birth-weight by ethnicity and smoking status of the mother.

5.4 Violin Plots

Violin plots are an alternative to box-plots.

plot(
  kfm,
  x = :sex, y = :dl_milk,
  Geom.violin,
  Guide.xlabel(""),
  Guide.ylabel("Breast milk intake (dl/day) "),
  Theme(default_color="Plum")
)

Figure 20: Comparison of breast-milk intake by sex.

plot(
  wcgs,
  x = :beh_pat, y = :sbp, color = :beh_pat, xgroup=:chd,
  Geom.subplot_grid(Geom.violin),
  Guide.xlabel("Coronary Heart Disease status"),
  Guide.ylabel("SBP (mmHg)"),
  Guide.colorkey("Behaviour \n Pattern")
)

Figure 21: Comparison of systolic blood pressure (SBP) by coronary heart disease event and behaviour pattern group.

5.5 Bar Charts

Bar charts or bar plots are used, most of the time, to display information from a table graphically. Therefore, for publications are not the best option as a table uses less space and contains better information.

5.5.1 Frequency

For categorical variables, we can plot the counts with frequency.

@combine(groupby(birth, [:race, :low]), :count = @nrow)

6×3 DataFrame

Row	race	low	count
	Cat…	Cat…	Int64
1	White	No	73
2	White	Yes	23
3	African-American	No	15
4	African-American	Yes	11
5	Other	No	42
6	Other	Yes	25

plot(
  @combine(groupby(birth, [:race, :low]), :count = @nrow),
  x = :race, y = :count, color = :low,
  Geom.bar(position=:dodge),
  Guide.xlabel(""),
  Guide.colorkey(title="Low birth weight?", pos=[0.4w, -0.3h])
)

Figure 22: Comparison of low-birth babies by ethnicity.

5.5.2 Expectation

For continuous variables, we can plot the expected (mean) values, with expectation.

combine(groupby(kfm, :sex), :dl_milk => cis => AsTable)

2×4 DataFrame

Row	sex	outcome	lower	upper
	Cat…	Float64	Float64	Float64
1	Boy	7.9524	7.36768	8.53712
2	Girl	7.0564	6.49906	7.61374

plot(
    combine(groupby(kfm, :sex), :dl_milk => cis => AsTable),
    x = :sex, y = :outcome, ymin = :lower, ymax = :upper,
    Scale.x_discrete(levels=["Boy", "Girl"]),
    Geom.errorbar, Stat.dodge, Geom.bar(position = :dodge),
    Guide.xlabel(""),
    Guide.ylabel("Breast-milk intake (dl/day)")
)

Figure 23: Comparison of daily breast-milk intake by sex.

combine(groupby(birth, [:smoke, :Race]), :bwt => cis => AsTable)

4×5 DataFrame

Row	smoke	Race	outcome	lower	upper
	Cat…	Cat…	Float64	Float64	Float64
1	Non-smoker	White	3428.75	3218.93	3638.57
2	Non-smoker	Non-white	2824.51	2664.85	2984.17
3	Smoker	White	2826.85	2656.57	2997.12
4	Smoker	Non-white	2642.09	2336.65	2947.54

plot(
    combine(groupby(birth, [:smoke, :Race]), :bwt => cis => AsTable),
    x = :Race, y = :outcome, color = :smoke, 
    ymin = :lower, ymax = :upper,
    Scale.x_discrete(levels=["White", "Non-white"]),
    Geom.errorbar, Stat.dodge, 
    Geom.bar(position = :dodge),
    Guide.xlabel(""),
    Guide.ylabel("Birth weight (g)"),
    Guide.colorkey(title="", pos=[0.8w, -0.4h])
  )

Figure 24: Comparison of birth weight by ethnicity and smoking status.